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Introduction 


This  report  describes  an  initial  construction  of  a  general  framework  for  numerical 
simulation  of  the  various  possible  types  of  scenarios  that  could  possibly  occur  for  the 
detection  and  remote  activation  of  improvised  explosive  devices  (IEDs)  by  excitation  of 
incident  electromagnetic  waves.  This  general  framework  consists  of  a  set  of  component 
models,  each  of  whose  structure  permits  the  output  of  given  types  of  information.  The 
model  representing  the  central  component  of  this  framework,  to  which  the  outputs  of  all 
the  other  component  models  are  inputs,  is  that  of  an  S-matrix  representation  of  a 
multilayered  composite  material  system,  where  each  layer  of  the  system  is  characterized 
by  an  average  thickness  and  effective  electric  permittivity  function  [1],  The  outputs  of  this 
primary  component  are  the  reflectivity  and  transmissivity  as  a  function  of  frequency  and 
incident  angle  of  the  incident  electromagnetic  wave.  The  other  component  models,  whose 
outputs  are  input  to  the  S-matrix  model,  are  response  spectra  calculated  using  density 
functional  theory  (DFT)  [2-4]  and  related  methodologies,  parameterized  analytic  function 
representations  of  the  electric  permittivity  as  a  function  of  frequency  obtained  by  fitting 
experimentally  measured  spectra,  and  effective  permittivity  functions  whose  construction 
is  based  on  effective  medium  theory  (EMT)  and  roughness  models.  We  review  those 
physical  theories  establishing  the  foundation  of  the  component  models  and  a  prototype 
simulation  that  considers  response  characteristics  for  THz  excitation.  We  include  an  initial 
version  of  a  computer  program  for  calculation  of  reflectivity  and  transmissivity  functions 
using  the  S-matrix  formulation.  Aspects  of  this  specific  software  implementation  are 
discussed.  In  addition,  we  describe  a  procedure  for  calculating  response  spectra  using  DFT 
for  use  as  input  to  the  S-matrix  model.  For  this  purpose  we  have  adopted  the  DFT 
software  NRFMOF. 

It  is  significant  to  note  that  the  numerical- simulation  framework  to  be  presented  is 
structured  for  two  major  purposes,  which  are  complimentary.  One  purpose,  which  relates 
directly  to  practical  application,  is  simulation  of  various  possible  scenarios  for  detection  of 
IEDs  corresponding  to  the  presence  of  various  types  of  intermediate  material  layers 
between  explosive  and  detector.  The  other  purpose,  which  relates  indirectly  to  practical 
application,  but  is  yet  extremely  important  for  the  interpretation  and  design  of  detection 
strategies,  is  the  quantitative  analysis  of  absolute  bounds,  or  rather,  the  inherent  limitation 
on  levels  of  detection  associated  with  various  types  of  detection  strategies.  With  respect  to 
the  purpose  of  examining  inherent  limitations  on  IED  detection,  the  dominant  features  of 
response  spectra  that  are  calculated  using  DFT  provide  a  foundation  for  establishing  what 
level  of  detection  is  achievable  in  the  absence  of  instrumental  and  environmental  factors 
associated  with  detection.  Accordingly,  the  simulation  framework  presented  here 
considers  a  specific  application  of  DFT.  For  any  given  energetic  material  and  frequency 
range  of  the  incident  electromagnetic  wave,  the  output  of  the  DFT  model  component  is  a 
set  of  response  signatures  that  are  each  characterized  by  an  excitation  frequency, 
magnitude  and  width.  These  response  signatures  must  then  be  used  to  construct 
permittivity  functions,  which  represent  the  form  of  input  to  the  S-matrix  component  of  the 
simulation  framework. 

A  significant  aspect  of  the  numerical- simulation  framework  presented  is  that  it 
adopts  the  perspective  of  computational  physics,  according  to  which  a  numerical 
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Figure  1.  General  framework  for  numerical  simulation  of  IED  detection  and  remote 
activation  scenarios. 


simulation  represents  another  source  of  “experimental”  data.  This  perspective  is 
significant  in  that  a  general  procedure  may  be  developed  for  construction  of  permittivity 
functions  that  uses  both  DFT  calculations  as  well  as  experimental  measurements.  That  is 
to  say,  for  the  purpose  of  simulating  many  electromagnetic  response  characteristics  of 
energetic  materials,  DFT  and  associated  methodologies  such  as  tight  binding  (TB) 
methods,  are  sufficiently  mature  for  the  purpose  of  generating  data  that  complements 
experimental  measurements  rather  than  simply  providing  verification. 


General  Simulation  Framework 

A  schematic  representation  of  the  general  framework  for  numerical  simulation  of 
IED  response  is  shown  in  Fig.  1.  It  should  be  emphasized  that  this  represents  an  initial 
construction  and  that  the  general  framework  shown  in  Fig.  1  is  subject  to  subsequent 
refinement  and  modifications  with  respect  to  the  paths  of  input  and  output  from  the 
different  model  components  comprising  the  framework.  Referring  to  Fig.  1,  it  should  be 
noted  that  the  primary  input  to  the  model  system  is  the  set  of  permittivity  functions  that 
are  associated  with  the  different  layers  of  material  making  the  system. 


Description  Of  Component  Models 

S-Matrix  Representation  of  Layered  Composite  System 

The  central  component  of  the  general  simulation  framework,  to  which  the  outputs 
of  all  the  other  component  models  are  inputs,  is  that  of  an  S-matrix  representation  of  a 
multilayered  composite  material  system,  where  each  layer  of  the  system  is  characterized 
by  an  average  thickness  and  effective  electric  permittivity  function  (see  [5]).  The  outputs 
of  this  central  component  are  the  reflectivity  and  transmissivity  as  a  function  of  frequency, 
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incident  angle  and  polarization  of  the  incident  electromagnetic  wave.  The  formulation  of 
the  S-matrix  representation  is  defined  by  the  following  system  of  equations. 

The  reflectivity  R  and  transmissivity  T  functions  are  given  by 

y  y  y  _  y  y 

R  =  -^  and  T=  11  22  12  21  ,  (1) 


respectively,  where  the  S-matrix  elements  Sl7  are  define  by  the  matrix  relation 


[S]  = 


Su  $ 


12 


V^21  S22,  J=l 


rn 

=ni“d 


where  [M j\  =  .  The  matrix  [/„,,]  is  defined  by  the  matrix  relation 


Lab 


\rab  1  / 


where 


r  _  £/A  ~£gSb 

ab  £bSa+£aSb 


and  tab  = 


2  S„ 


Sa+Sh 


for  a  p-polarized  incident  wave,  and 


rab  =  ^71T  and  tab  = 


2  ekS„ 


sa+sb 


£bSa  +£aSb 


for  an  s-polarized  incident  wave,  where 


2  j.\1/2 


sa  =(£a  -£0sin-(p)  -  and  Sb  =  {eb  -  e0  sin  (p) 


(2) 


(3) 

(4) 

(5) 

(6) 


Here,  £0  is  the  permittivity  function  of  the  transparent  ambient  layer,  £a  and  Eh  are  the 
permittivity  functions  for  layers  “a”  and  respectively.  The  matrix  [Ly]  is  defined  by 
the  matrix  relation 


where 


Ilj]  = 


‘xT 

0 


l,Xj 


1/2 


(7) 


f  dj  ) 

1 

[a] 

Xj  —  exp 

-2  iti 

(D^)} 

,  DjQ)  =  - 

=  (£.  -  £0  sin2  (p)v2 


(8) 
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and  Ej  is  the  permittivity  functions  for  layer  The  layer  indexing  used  in  Eqs.(l) 
through  (8)  is  defined  with  reference  to  Fig.  2.  A  computer  program  for  numerical 
implementation  of  Eqs.(l)  through  (8)  is  given  in  Appendix  1. 


Figure  2.  Schematic  representation  of  layer  indexing  used  in  Eqs.(l)  through  (8). 


Dielectric  Permittivity  Functions 

The  set  of  permittivity  functions  that  are  associated  with  the  different  layers 
comprising  the  layered  composite  system  represent  the  primary  input  to  the  S-matrix 
model  component.  These  functions  are  to  be  constructed  in  principle  for  a  given  material 
according  to  a  “best  fit”  of  available  information  associated  with  the  electromagnetic 
response  of  that  material.  For  a  given  material  this  information  consists  of  data  obtained 
from  both  experimental  measurements,  e.g.,  reflectivity  and  absorption  measurements,  and 
numerical  simulations  based  on  basic  principles,  e.g.,  DFT  and  TB  calculations.  It  is 
significant  to  note  that  the  best  fit  to  the  electromagnetic  response  of  a  given  material  will 
depend  on  the  specific  response  signature  characteristics  of  that  material.  Accordingly, 
from  the  perspective  of  numerical  simulation,  a  best  fit  can  be  in  the  form  of  a  tabulated 
functional  dependence,  as  well  as  representations  using  analytical  functions. 

There  are  specific  materials  that  are  typically  present  in  the  ambient  environment 
associated  with  IED  detection  as  well  as  the  detection  of  other  types  of  materials,  e.g., 
water  and  water  vapor.  Accordingly,  the  electromagnetic  response  characteristics  of  these 
materials  have  received  a  considerable  analysis  by  many  groups  and  are  available.  It 
follows  that  the  permittivity  functions  of  these  materials  should  represent  a  permanent 
“data  base”  component  of  the  general  simulation  framework.  An  example  of  the 
measurement  of  absorption  coefficients  of  selected  explosives  that  are  covered  by 
different  types  of  materials  (plastic,  cotton  and  leather)  are  given  in  Ref.(6). 
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Density  Functional  Theory  and  Related  Methodologies 

The  application  of  density  functional  theory  (DFT)  and  related  methodologies  for 
the  determination  of  electromagnetic  response  characteristics  is  important  for  the  analysis 
of  parameter  sensitivity.  That  is  to  say,  many  characteristics  of  the  electromagnetic 
response  of  a  given  material  may  not  be  detectable,  or  in  general,  not  relevant  for 
detection.  Accordingly,  sensitivity  analyses  concerning  the  electromagnetic  response  of 
layered  composite  systems  can  adopt  the  results  of  simulations  using  DFT,  and  related 
methodologies,  to  provide  realistic  limits  on  detectability  that  are  independent  of  a  specific 
system  design  for  IED  detection.  In  addition,  analysis  of  parameter  sensitivity  based  on 
atomistic  response  characteristics  of  a  given  material,  obtained  by  DFT,  provide  for  an 
“optimal”  best  fit  of  experimental  measurements  for  the  construction  of  permittivity 
functions.  It  follows  that  within  the  context  of  parameter  sensitivity  analysis,  data  obtained 
by  means  of  DFT  represents  a  true  complement  to  data  that  has  been  obtained  by  means  of 
experimental  measurements. 

Experimental  Measurements 

The  dominant  amount  of  information  that  is  adopted  for  the  construction  of 
permittivity  functions  is  obtained  from  experimental  measurements  of  electromagnetic 
response  characteristics.  Some  major  issues  associated  with  these  constructions  are  that 
such  experimental  measurements  typically  involve  bulk  material  response  characteristics 
as  well  as  measurement  errors  due  to  sample  surface  preparations  and  artifacts  due  to 
ambient  environmental  influences.  These  issues  are  significant  in  that  the  permittivity 
functions  adopted  as  input  are  typically  assumed  as  being  associated  with  “pure”  materials 
as  well  as  representative  of  response  characteristics  on  a  small  scale  that  may  be  typical  of 
thin  film  type  layers.  As  in  the  case  of  response  characteristics  that  are  determined  via 
atomistic  calculations,  certain  response  features  associated  with  response  characteristics 
determined  by  experimental  measurement  may  not  be  significant  for  the  simulation  of  IED 
detection.  That  is  to  say,  certain  features  such  as  the  locations  and  amplitudes  of  response 
spectra  may  be  significant  for  inclusion  into  model  representations,  while  only  a 
reasonable  estimate  of  the  widths  may  be  necessary.  It  follows  that  sensitivity  analysis  for 
parameterizations  of  experimental  measurements  is  as  relevant  as  those  associated  with 
theoretical  predictions.  Such  analysis  is  another  application  goal  of  the  simulation 
framework  presented  here. 

Effective  Medium  Theory  (EMT),  Equivalent  Layers  and  Roughness  Models 

Consistent  with  the  goal  of  determining  absolute  limitations  on  detectability  of 
IEDs  by  means  of  electromagnetic  excitation  is  the  construction  of  models  of  material 
response  that  are  representative  of  a  general  class  of  materials  and  detection  scenarios,  in 
contrast  to  models  that  would  tend  to  be  associated  with  a  specific  experimental 
arrangement  in  the  laboratory.  Accordingly,  the  concepts  of  an  “effective  medium”  and 
“equivalent  layer”  are  significant  in  that  their  consideration  for  model  construction  can 
provide  quantitative  bounds  on  detectability  for  a  wide  range  of  detection  scenarios.  In 
particular,  these  concepts  can  provide  a  foundation  for  the  parametric  representation  of 
surface  roughness  and  inhomogeneities  on  various  spatial  scales  in  the  ambient 
environment. 
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The  formal  structure  of  all  continuum  effective  medium  theories  (EMTs)  are  based 
on  mixing  rules  that  are  functions  of  the  different  permittivities  making  up  the  composite 
material  (See  [7]  and  references  therein  for  further  discussion  of  EMTs). 

The  concept  of  an  equivalent  layer  is  based  on  the  fact  that  a  given  range  of 
different  types  of  layer  structures  can  have  the  same  average  response  characteristics. 
Accordingly,  a  layer  structure  that  is  within  this  range  can  be  represented  by  means  of  an 
equivalent  layer  whose  construction  does  not  require  consideration  of  many  details 
associated  with  its  composition.  A  particular  case  of  the  application  of  the  concept  of  an 
equivalent  layer  in  combination  with  EMT  is  the  construction  of  roughness  models  for  the 
representation  of  rough  surface  structure. 


Prototype  Analysis  (THz  Excitation) 

Presented  in  this  section  is  a  prototype  simulation  for  demonstrating  some  aspects 
of  the  relationship  between  the  various  model  components  that  comprise  the  general 
simulation  framework.  For  this  simulation  the  response  of  a  layer  of  /3-HMX  to  THz 
excitation  is  considered  [8-10]. 

Shown  in  Fig.  3  are  the  real  and  imaginary  parts  of  a  permittivity  function 
corresponding  to  the  electromagnetic  response  of  /3-HMX  to  excitation  within  the  THz 
range  of  frequencies.  This  permittivity  function  is  significant  in  that  it  has  been 
constructed  using  DFT  calculations  that  have  been  calibrated  with  reference  experimental 
measurements.  Accordingly,  the  approach  followed  for  construction  of  this  permittivity 
function  is  that  which  has  been  adopted  for  the  construction  of  permittivity  functions 
within  the  simulation  framework  presented,  i.e.,  a  best  fit  to  the  combination  of  both 
theoretical  calculations  and  experimental  data. 

Shown  in  Fig.  4  are  reflectivity  functions  corresponding  to  s  and  p-polarization  of 
the  incident  wave.  The  layered  system  consists  of  a  layer  of  /3-HMX  upon  a  gold 
substrate.  The  reflectivity  functions  shown  in  Fig.  4,  in  principle,  would  represent  the 
starting  point  for  any  study  concerning  absolute  bounds  on  the  detectability  of  /3-HMX 
under  different  environmental  conditions  (i.e.,  surface  layers  and  ambient  environment) 
and  detection  scenarios. 
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Figure  3.  Permittivity  function  of  /j-HMX  for  frequencies  within  THz  range. 


Figure  4.  Reflectivity  functions  for  a  layer  of  /3  -HMX  on  a  gold  substrate. 
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Analysis  Of  Permittivity  Functions  Using  DFT 


The  general  approach  of  constructing  permittivity  functions  according  to  the  best 
fit  of  available  data  for  given  material  corresponding  to  many  different  types  of 
experimental  measurements  is  not  unprecedented  and  has  been  typically  the  dominant 
approach,  e.g.,  the  permittivity  function  of  water.  The  general  simulation  framework 
presented  here  considers  an  extension  of  this  approach  in  that  calculations  of 
electromagnetic  response  based  on  DFT  and  associated  methodologies  are  also  adopted  as 
data  for  construction  of  permittivity  functions.  The  inclusion  of  this  type  of  information  is 
significant  for  accessing  what  spectral  response  features  at  the  molecular  level  are  actually 
detectable  with  respect  to  a  given  set  of  detection  parameters.  Accordingly,  permittivity 
functions  having  been  constructed  using  DFT  calculations  provide  a  quantitative 
correlation  between  macroscopic  material  response  and  molecular  structure.  Within  this 
context  it  is  not  important  that  the  permittivity  function  be  quantitatively  accurate  for  the 
purpose  of  being  adopted  as  input  for  system  simulation.  Rather,  it  is  important  that 
permittivity  function  be  qualitatively  accurate  in  terms  of  its  general  features  for  the 
purpose  of  sensitivity  analysis,  which  is  relevant  for  the  assessment  of  absolute 
detectability  of  different  types  of  molecular  structure  with  respect  to  a  given  set  of 
detection  parameters.  That  is  to  say,  permittivity  functions  that  have  been  determined 
using  DFT  can  provide  a  mechanistic  interpretation  of  material  response  to 
electromagnetic  excitation  that  could  establish  the  well  posedness  of  a  given  detection 
methodology  for  detection  of  specific  molecular  characteristics.  Within  the  context  of 
practical  application,  permittivity  functions  having  been  constructed  according  to  the  best 
fit  of  available  data  would  be  “correlated”  with  those  obtained  using  DFT  for  proper 
interpretation  of  permittivity-function  features.  Subsequent  to  establishment  of  good 
correlation  between  DFT  and  experiment,  DFT  calculations  can  be  adopted  as  constraints 
for  the  purpose  of  constructing  permittivity  functions,  whose  features  are  consistent  with 
molecular  level  response,  for  adjustment  relative  to  specific  sets  of  either  experimental 
data  or  additional  molecular  level  information. 


Figure  5.  General  procedure  for  construction  of  permittivity  functions  using  DFT 
calculations. 

The  construction  of  permittivity  functions  using  DFT  calculations  involves, 
however,  an  aspect  that  requires  serious  consideration.  This  aspect  concerns  the  fact  that  a 
specific  parametric  function  representation  must  be  adopted.  This  significant  aspect  of 
constructing  permittivity  functions  using  DFT,  and  related  methodologies,  is  shown 
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explicitly  in  Fig.  5.  Accordingly,  any  parametric  representation,  i.e.,  parameterization, 
adopted  for  permittivity-function  construction  must  be  physically  consistent  with  specific 
molecular  response  characteristics,  while  limiting  the  inclusion  of  feature  characteristics 
that  tend  to  mask  response  signatures  that  may  be  potentially  detectable. 

In  principle,  parameterizations  are  of  two  classes.  One  class  consists  of 
parameterizations  that  are  directly  related  to  molecular  response  characteristics.  This  class 
of  parameterizations  would  include  spectral  scaling  and  width  coefficients.  The  other  class 
consists  of  parameterizations  that  are  purely  phenomenological  and  are  structured  for 
optimal  and  convenient  best  fits  to  experimental  measurements. 

At  this  stage  it  is  instructive  to  present  a  prototype  calculation  demonstrating 
analysis,  e.g.,  interpretation,  of  permittivity-function  features  using  DFT  calculations. 
Consistent  with  the  prototype  simulation  presented  above,  a  permittivity  function  is 
constructed  using  DFT  calculations  for /3-HMX  response  to  THz  excitation.  Shown  in 
Fig.  6  is  a  general  description  of  the  geometry  of  the  /3-HMX  molecule  that  was  adopted 
for  calculation  of  a  permittivity  function  using  DFT.  That  is  to  say,  the  molecular  structure 
that  was  input  to  the  DFT  software  NRLMOL.  Shown  in  Figs.  7  and  8  are  absorption 
coefficients  corresponding  to  different  adjustable  spectral  scaling  and  width  parameters. 


Figure  6.  Molecular  structure  of  /3-HMX  used  for  DFT  calculations  of  spectral  response. 

Remark.  It  is  significant  to  note  that  with  respect  to  practical  application  the 
absorption  coefficient  a 


a  = 


4  n 


V2A 


-£r  + 


4 


e:+£7 


1/2 


(9) 
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where  A,  er  and  sj  are  the  wavelength  of  excitation,  and  the  real  and  imaginary  parts  of 
the  permittivity  function,  respectively,  provides  a  direct  relationship  between  a  calculated 
quantity  using  DFT  and  a  “conveniently  measurable”  quantity  a. 

Next,  we  consider  a  qualitative  example  of  examining  the  correlation  between 
DFT  calculated  permittivity  functions  and  experimental  measurements.  Referring  to  Fig. 
9,  which  shows  an  experimentally  determined  absorption  coefficient  for  j>  -  HMX  (see 
Ref.(7)),  we  note  good  correlation  between  the  permittivity  functions  (in  terms  of  their  a 
representation)  obtained  by  DFT  (using  NRLMOL)  and  experiment.  Most  importantly,  the 
level  of  correlation  is  sufficient  to  establish  a  “proof  of  concept”  that  DFT  calculations 
provide  a  quantitative  initial  estimate  of  molecular  response  to  electromagnetic  excitation 
for  subsequent  parameterization  [11-13]. 


Figure  7.  Absorption  coefficient  for  /J  -HMX  calculated  by  DFT  for  THz  range  of 
frequencies  corresponding  to  adjustable  parameters  y  =  3  cm  1  and  N  =  2  cm  2. 

As  indicated  previously,  the  parameterizations  applied  to  DFT  calculations  will  in 
general  consist  of  two  classes  of  parameterizations,  i.e.,  one  consistent  with  basic  theory 
and  the  other  consistent  with  optimal  and  convenient  best  fitting  of  experimental 
measurements.  Accordingly,  one  class  of  parameterization  defines  a  problem  requiring 
further  analysis  in  terms  of  basic  theory  [8,10,14,15],  while  the  other  class  defines  a 
problem  requiring  analysis  in  terms  of  inverse-problem  and  parameter-optimization 
methodologies  [16]. 
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Figure  8.  Absorption  coefficient  for  /8-HMX  calculated  by  DFT  for  THz  range  of 
frequencies  corresponding  to  adjustable  parameters  y  =  5  cm  1  and  N  =  3  cm  2. 


Figure  9.  Experimentally  determined  absorption  coefficient  for  [>  -  HMX  for  THz  range 
of  frequencies  (see  [8]). 
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Conclusion 


The  development  of  IED  detection  methodologies  requires  the  consideration  of 
two  major  aspects  of  detection.  These  are  the  detectability  of  IEDs  under  different  types  of 
environmental  conditions  and  detection  scenarios;  and  the  absolute  detectability  of  the 
different  types  of  response  characteristics  of  energetic  materials  due  to  electromagnetic 
wave  excitation.  Accordingly,  within  the  context  of  practical  application  of  IED  detection 
methodologies,  it  remains  necessary  to  establish  correlation  with  response  properties  on 
the  molecular  level.  It  is  therefore  necessary  to  construct  two  types  of  permittivity 
functions.  One  type,  whose  purpose  is  the  simulation  of  detection  scenarios,  represents  the 
best  fit  to  available  data,  which  could  include  both  experimental  measurements  and 
calculations  based  on  theory.  The  other  type,  obtained  using  DFT,  is  that  of  a  reasonably 
optimal  parametric  representation  of  molecular  level  response  characteristics,  providing 
interpretation  of  permittivity-function  features  at  the  molecular  level,  whose  purpose  is 
that  of  initial  constraints  for  subsequent  adjustment  relative  to  specific  sets  of  either 
experimental  data  or  additional  molecular  level  information.  It  follows  that  the 
establishment  of  a  general  constrained  parameterization  of  IED  response  based  on  both 
theory  and  experiment,  combined  with  quantitative  sensitivity  analyses  of  this 
parameterization,  will  provide  an  assessment  of  the  general  detectability  of  IEDs 
independent  of  specific  detection  scenarios. 
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Disclaimer 

This  report  does  not  imply  any  form  of  warranty  that  the  S-Marix  code  does  not  contain 
errors  or  that  it  is  sufficient  for  any  specific  application.  The  S-Matrix  code  should  not  be 
relied  on  for  solving  problems  whose  incorrect  solution  could  result  in  damages. 
Accordingly,  the  authors  of  this  S-Matrix  code  and  this  report  disclaim  all  liability  for 
direct  or  consequential  damages  resulting  from  the  use  of  the  S-Matrix  code. 
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Appendix  1 


Computer  Program  For  Calculation  Of  S-Matrix 

Presented  in  this  section  is  an  initial  version  of  a  computer  program,  i.e.,  the  “S- 
Matrix  code”  (in  Fortran  77)  for  calculation  of  the  reflectivity  and  transmissivity  function 
using  the  S-Matrix  representation  of  a  multilayered  composite  material  system. 


PROGRAM  MLAYERSP 
C 

C  version  21 

C  last  time  modified:  03/23/2010 

C 

implicit  none 

C  parameters  and  constants 

integer  ML,LNR,NMP  !  max  number  of  layers,  resonances,  mesh  points 
parameter  (ML  =  20,  LNR  =  100,  NMP  =  10000) 

real*8  c,pi  !  speed  of  light  and  pi 

parameter  (c  =  2 . 99792458d8,  pi=3 . 14159265358979323846d0) 


complex*16  SP (2, 2) , SS  (2, 2)  !  scattering  matrix 

real*8  theta  (NMP  ),  angle  (NMP  )  !  incident  angle 

real*8  ang0,angl  !  and  its  range 

real*8  wnumO , wnuml  !  range  of  wave  numbers 

real*8  omega  (NMP ),  lamda  (NMP )  !  angular  frequency  and  wave  length 

real  *8  wnum  (NMP  )  ,  awnum  (NMP  )  !  wavenumber  and  angular  wavenumber 

COmplex*16  E0 (NMP) , E (ML, NMP) , EREAD (NMP) 

real*8  D (ML)  !  thickness  of  each  layer 

COmplex*16  REFP , REFS, TRANP , TRANS, DETSS, DETSP 

real*8  ref lb, ref2b 

real*8  ref 1 (NMP ), ref 2 (NMP ) 

integer  IUNITS  !  input  units 

integer  ISCAN  !  type  of  scan 

integer  NLAYER,  NFACE,  NREGION  !  number  of  layers,  faces,  regions 
integer  fol  !  first  opaque  layer 

integer  NMESHA, NMESHF, NBINA, NBINF  !  number  of  angles,  freq,  binsize 
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integer  NBINA2 , NBINF2  !  half  of  binsize 


integer  IL, I, J, IB, JB, ILS  !  counters 

character*79  layerfl  !  layer  data  file 
integer  ecode  !  errorcode 


open  ( 1 ,  f i le=  ' MLAYER .  INP  '  )  !  open  input  and  output  files 

open (2, file= 'MLAYER. OUT '  ) 
open  (3, f ile=' MLAYER. AVG' ) 

read ( 1 , * )  IUNITS,  ISCAN  !  input  units  and  type  of  scan 
read(l,*)  wnumO , wnuml , NMESHF, NBINF  !  wave  numbers 
NBINF  =  2*INT (  (NBINF  +  1) /  2) -1 
NBINF2=INT (NBINF/ 2 ) 
do  J=1 , NMESHF +NBINF-1 

wnum(J)  =  wnumO  +  (float (J-l) /float (NMESHF) )* (wnuml-wnumO ) 

lamda ( J) =1 . OdO/wnum ( J)  !  wave  number 

awnum  ( J)  =2  .  OdO *pi *wnum  ( J)  !  angular  wave  number 

if  (IUNITS  .eq.  0)  then 

omega (J)  =  1 . 0d2*awnum ( J) *c 
else 

omega (J)  =  awnum (J)*c 
end  if 
end  do 

read ( 1 , * )  angO,  angl,  NMESHA, NBINA  !  angles 
NBINA=2  * INT (  (NBINA+1) / 2) -1 
NBINA2=INT (NBINA/ 2 ) 
do  1=1 , NMESHA+NBINA-1 

theta  (I)  =  angO  +  ( float ( 1-1 ) /float (NMESHA) )* (angl-angO ) 
theta(I)  =  (theta ( I ) / 1 . 8d2 ) *pi 
angle(I)  =  (theta (I) /pi) *1 . 8d2 
end  do 

read(l,*)  NLAYER  !  number  of  layers 
NFACE  =  NLAYER  +  1  !  number  of  interfaces 

NREGION  =  NLAYER  +  2  !  number  of  regions 

read(l,*)  !  blank  line 
read(l,  '  (a79)  ')  layerfl 

call  READ LAYER ( layer f 1 , NMESHF , wnum, EREAD , ecode )  !  ambient  layer 

if  (ecode  .eq.  0)  then 
do  J=l, NMESHF 

E0 ( J)  =  EREAD ( J) 
end  do 
else 

print  *,  'INPUT  ERROR  =  1 , ecode 
print  *,  'CANNOT  READ  AMBIENT  LAYER' 
stop 
end  if 

do  IL=1, NFACE  !  read  layers 
print  *,  ' LAYER= ' , I L 

read(l,*)  !  blank  line 

read(l,  '  (a79)  ')  layerfl  !  layer  data  file  name 
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read(l,*)  D  ( I L )  !  thickness 

call  READ LAYER (layer f 1, NMESHF, wnum, EREAD, ecode) 
if  (ecode  .eq.  0)  then 
do  J=l, NMESHF 

E ( IL, J)  =  EREAD ( J) 
end  do 
else 

print  *,  'INPUT  ERROR  =  ecode 
print  *,  'CANNOT  READ  LAYER  =  ' , IL 

stop 
end  if 

if  (D(IL)  .It.  4.0d0/wnuml)  then 

print  *,  'WARNING!  thickness  is  too  small' 
end  if 

end  do  !  end  of  layers 

if  (ISCAN  .eq.  0)  then  !  angle  scan 

write (2,*)  '  ANGLE  Rp  Rs  ' 

do  I  =  1 , NMESHA+NBINA-1 
ref lb=0 . OdO 
ref 2b=0 . OdO 
do  J  =  1 , NBINF 

ILS=0  !  set  control  layer  number  ILS  and  fol  to  0 
fol=0 

do  IL  =  1 , NFACE 

if  (ILS  .It.  NFACE)  then  !  call  layers  if  JS  <  NFACE 
f ol=IL 

if  (IL  .eq.  1)  then  !  initial  call  for  ambient  layer 
call  LAYER ( IL, NLAYER, E0,E0,E(IL,J),  theta ( I)  , 

&  lamda ( J) , 0 . OdO, SP, SS, ILS) 

else 

call  LAYER ( IL, NLAYER, E0,E(IL-1,J),E(IL,J), theta ( I) , 
&  lamda (J)  ,  D (IL-1) , SP, SS, ILS) 

end  if  !  IL=1 
end  if  !  ILS<NFACE 
end  do  !  end  of  layers 
REFP  =  -SP  (1, 2) /SP  (1,  1) 

DETSP=SP (1, 1) *SP (2,2) -SP (1,2)*SP(2,1) 

TRANP  =  DETSP/SP  (1,  1) 

REFS  =  -SS  (2,  1) /SS  (1,1) 

DETSS=SS (1, 1) *SS (2, 2) -SS (1,2)*SS(2,1) 

TRANS  =  DETSS/SS  (1,  1) 
ref lb=ref lb+REFP*con jg (REFP) 
ref 2b=ref lb+REFS*con jg (REFS) 
end  do  !  end  of  frequency  bin 
ref 1 ( I ) =re fib /NBINF 
ref 2 ( I ) =ref 2b /NBINF 

write(2,*)  angle (I),  refl(I),  ref2(I) 
end  do  !  end  of  angle  scan 

write (3,*)  '  ANGLE  Rp  Rs  ' 

do  I  =  1 , NMESHA 
ref lb=0 . OdO 
ref 2b=0 . OdO 
do  IB  =  1 , NBINA 

ref lb=ref lb+ref 1 ( I-NBINA2+IB-1 ) 
ref 2b=ref 2b+ref 2 ( I-NBINA2+IB-1 ) 
end  do 
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ref lb=ref lb/NBINA 
ref 2b=ref 2b/NBINA 

write  (3,*)  angle (J),  reflb,  ref2b 
end  do 

else  !  frequency  scan 

write (2,*)  '  wnum  frequency  Rp  Rs  ' 

do  J  =  1 , NMESHF+NBINF-1 
ref lb=0 . OdO 
ref 2b=0 . OdO 
do  I  =  1 , NBINA 

ILS=0  !  set  control  layer  number  JS  and  fol  to  0 
f  ol=0 

do  IL  =  1 , NFACE 

if  (ILS  .It.  NFACE)  then  !  call  layers  if  ILS  <  NFACE 
f ol=IL 

if  (IL  .eq.  1)  then 

call  LAYER (IL, NLAYER, E0,E0,E(IL,J),  theta ( I)  , 

&  lamda ( J) , 0 . OdO, SP, SS, ILS) 

else 

call  LAYER (I L, NLAYER, E0,E(IL-1,J),E(IL,J), theta (I) , 
&  lamda (J)  ,  D (IL-1) , SP, SS, ILS) 

end  if  !  IL=1 
end  if  !  JS<NFACE 
end  do  !  end  of  layers 
REFP  =  SP (1, 2) /SP  (1,1) 

DETSP=SP (1, 1) *SP (2,2) -SP (1,2)*SP(2,1) 

TRANP  =  DETSP/SP  (1,  1) 

REFS  =  SS (1,2) /SS  (1,1) 

DETSS  =  SS (1, 1) *SS (2, 2) -SS  (1,2)*SS(2,1) 

TRANS  =  DETSS/SS  (1,  1) 
ref lb=ref lb+REFP*con jg (REFP) 
ref 2b=ref 2b+REFS*con jg (REFS) 
end  do  !  end  of  angle  bin 
refl ( J) =ref lb/NBINA 
ref2 ( J) =ref2b/NBINA 

write (2,*)  wnum(J),  omega (J),  refl (J) ,  ref2(J),  fol 
end  do  !  end  of  frequency  scan 

write (3,*)  '  wnum  frequency  Rp  Rs  ' 

do  J  =  1 , NMESHF 
ref lb=0 . OdO 
ref 2b=0 . OdO 
do  JB  =  1 , NBINF 

ref lb=ref lb+ref 1 ( J-NBINF2+ JB-1 ) 
ref 2b=ref 2b+ref 2 ( J-NBINF2+ JB-1 ) 
end  do 

ref lb=ref lb/NBINF 
ref 2b=ref 2b /NBINF 

write (3,*)  wnum(J),  omega (J),  reflb,  ref2b 
end  do 

end  if  !  angle  or  frequency  scan 

close  (1)  !  close  input/output 

close  (2 ) 
close  ( 3 ) 

stop 
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end  !  end  of  program 


C - SUBROUTINES - 

subroutine  NIM2ALPHA (nim, k, alpha) 

C  converts  absorption  coefficient  to  imaginary  part  of  refraction  index 

implicit  none 

real*8  nim,  alpha, k  !  im  part  of  ref  index,  abs  coef,  ang  wavenumber 

alpha  =  (2 . OdO*k) *nim 

return 

end 


subroutine  ALPHA2NIM (alpha, k, nim) 

C  converts  absorption  coefficient  to  imaginary  part  of  refraction  index 

implicit  none 

real*8  nim,  alpha,  k  !  im  part  of  of  ref  index,  abs  coef,  ang  wavenumber 

nim  =  alpha/ (2.0d0*k) 

return 

end 


subroutine  N2E (nre, nim, ere, eim) 

C  converts  index  of  refraction  to  dielectric  constant 

implicit  none 
real*8  nre, nim, ere, eim 
ere  =  nre**2  -  nim**2 
eim  =  2 . OdO*nre*nim 
return 
end 


subroutine  E2N (nre, nim, ere, eim) 

C  converts  index  of  refraction  to  dielectric  constant 

implicit  none 
real*8  sqrt 
intrinsic  sqrt 
real*8  nre, nim, ere, eim 

nre  =  sqrt (ere  +  sqrt (ere**2  +  eim**2 ) ) /sqrt (2 . OdO ) 
nim  =  sqrt (-ere  +  sqrt (ere**2  +  eim**2 ) ) /sqrt (2 . OdO ) 
return 
end 


subroutine  SUMOFRES (ere, eim, einf , vol, wnum, 

&  wnumN, sxN, syN, szN, gN, NRES) 

C  calculates  dielectric  constants 

implicit  none 
real*8  sqrt 
intrinsic  sqrt 
real*8  pi 

parameter  (pi=3 .1415 92 6535897932384 6d0) 
real *8  ere, eim, einf, vol , wnum, wnumCM 
real *8  wnumN (*) , sxN (*) , syN (*) , szN (*) , gN (*) 
real*8  stot,pref 
integer  NRES,  I 
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pref = (4  *pi/vol ) 
ere=einf 
eim=0 . OdO 
do  1  =  1 , NRES 

stot= ( sxN (I) +syN (I) +szN (I) ) /3 . OdO 

ere  =  ere  +  pref*stot* (wnumN (I) **2-wnum**2) / 

&  ( (wnum**2-wnumN (I) **2 ) **2+gN (I) **2*wnum**2 ) 

eim  =  eim  +  pref *stot* (wnum*gN ( I ) ) / 

&  ( (wnum**2-wnumN (I)**2)**2+gN(I) **2*wnum**2 ) 

end  do 
return 
end 

subroutine  LAYER (JL, JN, EO, elN, eOUT, theta, lamda,  D,  SP,  SS,  JS) 

C -  calculates  reflection  and  transmission 

C  of  em  wave  through  a  single  layer  interface 

C  and  accumulative  reflection  and  transmission 

C  layer  k-1 

C  -  interface 

C  layer  k 

- input /output  parameters 

JL  -  layer  number  (=  k) 

JN  -  number  of  layers  (NLAYER) 

EO  -  dielectric  constant  of  vacuum  (abandoned  in  ver.  6) 

elN  -  dielectric  constant  of  layer  k-1 

eOUT  -  dielectric  constant  of  layer  k 

theta  -  incident  angle 

lamda  -  wavelength 

D  -  thickness  of  layer  k 

SP,SS  -  scattering  matrices 

JS  -  flag  (=JL  for  transparent  layers)  set  to  NFACE  for  opaque  layer 
- internal  paramters 

RP,RS,TP,TS  -  complex  reflection  and  transmission  coefficients 

MS,MP  -  layer  k  matrices 

IP, IS  -  interface  matrices 

NFACE  -  NLAYER+ 1 

OD  -  optical  density  of  layer 

C  JT  -  flag  (=1  for  transparent  layer)  set  to  0  for  opaque  layer 

implicit  none 
real*8  pi 

parameter  (pi=3 .1415 92 6535897932384 6d0) 

COMPLEX* 1 6  eO , elN, eOUT 
real*8  theta, lamda,  D 

complex*16  RP,RS,TP,TS  !  complex  ref  and  trans  coefficients 

complex*16  MS (2, 2),  MP(2,2)  !  layer  matrices 

complex*16  L(2,2) 

complex*16  IP (2, 2) , IS  (2, 2) 
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complex*16  SP(2,2),  SS(2,2)  !  scattering  matrices 


C  complex*16  SA,  SB,  EA,  EB,  DB,  DJ,  XJ,  YJ 

COmplex*16  SA,  SB,  EA,  EB,  DA,  DJ,  XJ,  YJ 
complex*16  tSP  (2 , 2 ) , tSS  (2 , 2 ) 
complex*16  iunit,  A(2,2) 

complex*16  OD  !  optical  density 

integer  JL,  JN, JS, NFACE, JT 


C 


C 


c 


c 


c 


c 


c 


real*8  exp,  sqrt ,  sin,  dble  !  external  functions 
intrinsic  exp, sqrt, sin, dble 
C0MPLEX*16  cdsqrt , cdexp, dcmplx 
intrinsic  cdsqrt , cdexp, dcmplx 


initialize 

tSP/tSS  with  input  matrix 

tSP(l,l)  = 

SP  (1, 1) 

t SP (2,2)  = 

SP  (2,2) 

tSP (1,2)  = 

SP  (1, 2) 

tSP (2,1)  = 

SP  (2,1) 

tSS(l,l)  = 

SS  (1, 1) 

t SS ( 2 , 2 )  = 

SS  (2, 2) 

tSS (1,2)  = 

SS  (1, 2) 

t SS ( 2 , 1 )  = 

SS (2,  1) 

set  JS  and  JT  flags  to  'transparent'  values 

JS= JL 

JT=1 

initilaize  internal  parameters 
NFACE=  JN  +  1 

iunit  =  dcmplx ( 0 . OdO , 1 . OdO ) 

EA  =  elN 
EB  =  eOUT 

SA  =  cdsqrt (EA  -  EO * ( s in (thet a ) ) * *2 ) 

SB  =  cdsqrt (EB  -  EO *  ( s in  (thet a ) ) * *2 ) 


RP  =  (EB*SA  -  EA*SB) / (EB*SA  +  EA*SB) 
RS  =  (SA  -  SB) / (SA  +  SB) 

TP  =  2 . *EB*SA/ (EB*SA  +  EA*SB) 

TS  =  2 . *SA/ (SA  +  SB) 


Calculate 

IP11(1),  IP12  ( 1 ) ,  IP21  ( 1 ) , 

IP(1,1)  = 

1  .  /TP 

IP ( 1 / 2 )  = 

RP/TP 

IP (2, 1)  = 

IP  (1,2) 

IP (2,2)  = 

IP  (1, 1) 

Calculate 

ISll(l),  IS12  ( 1 ) ,  IS21  (1) , 

IS(1,1)  = 

1  .  /TS 

IS ( 1 , 2 )  = 

RS/TS 

IS (2, 1)  = 

IS (1, 2) 

IS ( 2 , 2 )  = 

IS (1, 1) 

DA  =  dcmplx (D, 0 . OdO ) 

DJ  =  ( lamda/2 . ) * (1 . /SA) 

OD  =  -2 . *iunit*pi* (DA/DJ) 


IP22  (1) 


IS22  (1) 


20 


c 

if  (DBLE(OD)  .gt.  1.0d2)  then  !  TEST  OPTICAL  DENSITY 
JS=NFACE 
JT=0 
end  if 

if  (JL  .eq.  1)  then  !  generate  initial  SP  and  SS 
tSP(l,l)  =  dcmplx ( 1 . OdO , 0 . OdO ) 
tSP(2,2)  =  dcmplx ( 1 . OdO , 0 . OdO ) *JT 
tSP(l,2)  =  dcmplx ( 0 . OdO , 0 . OdO ) 
tSP(2,l)  =  dcmplx ( 0 . OdO , 0 . OdO ) 
tSS(l,l)  =  dcmplx ( 1 . OdO , 0 . OdO ) 
tSS(2,2)  =  dcmplx ( 1 . OdO , 0 . OdO ) 
tSS(l,2)  =  dcmplx ( 0 . OdO , 0 . OdO ) 
tSS(2,l)  =  dcmplx ( 0 . OdO , 0 . OdO ) 
end  if 


if  (JT  .eq.  1)  then  !  'optically  thin'  layer 
C  Calculate  Lll(l),  L12(l),  L21(l),  L22(l) 

YJ  =  cdexp (0 . 5dO*OD) 

C  Calculate  propagation  matrix 

L  (1, 1)  =  YJ 

L  (2, 2)  =  1 . /YJ 
L(l,2)  =  dcmplx ( 0 . OdO , 0 . OdO ) 

L(2,l)  =  dcmplx ( 0 . OdO , 0 . OdO ) 

C  Calculate  MP11(1),  MP12(1),  MP21(1),  MP22(1) 

MP  (1,1)  =  IP  (1,  1)  *L  (1,  1)  +  IP  (1,2)  *L  (2,  1) 

MP  (2,1)  =  IP  (2,  1)  *L  (1,  1)  +  IP  (2, 2)  *L  (2,  1) 

MP  (1,2)  =  IP  (1, 1)  *L  (1, 2)  +  IP  (1, 2)  *L  (2, 2) 

MP  (2,2)  =  IP  (2, 1) *L  (1,2)  +  IP  (2, 2) *L  (2, 2) 

C  Calculate  MS11(1),  MS12(1),  MS21(1),  MS22(1) 

MS  (1,1)  =  IS (1, 1) *L (1, 1)  +  IS (1,2) *L(2, 1) 

MS  (2,1)  =  IS  (2, 1) *L  (1, 1)  +  IS  (2,2) *L(2, 1) 

MS  (1,2)  =  IS  (1, 1) *L  (1, 2)  +  IS  (1, 2) *L  (2, 2) 

MS  (2,2)  =  IS  (2, 1) *L  (1,2)  +  IS (2 , 2 ) *L  (2 , 2 ) 

C  assign  MP  matrix  to  em  matrix 

A  ( 1 , 1 )  =  MP (1, 1) *tSP (1, 1) +MP (1,2) *tSP  (2, 1) 

A  ( 2 , 1 )  =  MP (2, 1) *tSP (1, 1) +MP (2,2) *tSP (2, 1) 

A  ( 1 , 2 )  =  MP (1, 1) *tSP (1, 2) +MP (1,2) *tSP (2,2) 

A  ( 2 , 2 )  =  MP (2, 1) *tSP (1, 2) +MP (2,2) *tSP (2,2) 

tSP  ( 1 , 1 )  =  A ( 1 , 1 ) 

tSP  (1,2)  =  A ( 1 , 2 ) 

tSP  (2,1)  =  A(2,l) 

tSP  (2,2)  =  A (2, 2) 

C  assign  MS  matrix  to  em  matrix 

A  ( 1 , 1 )  =  MS (1, 1) *tSS (1, 1) +MS (1,2) *tSS (2, 1) 

A  ( 1 , 2 )  =  MS (1, 1) *tSS (1, 2) +MS (1,2) *tSS (2,2) 

A  ( 2 , 1 )  =  MS (2, 1) *tSS (1, 1) +MS (2,2) *tSS (2, 1) 

A  ( 2 , 2 )  =  MS (2, 1) *tSS (1, 2) +MS (2,2) *tSS  (2,2) 
tSS  (1,1)  =  A(l,  1) 
tSS  (1,2)  =  A ( 1 , 2 ) 
tSS  (2,1)  =  A(2, 1) 
tSS  (2,2)  =  A (2, 2) 
end  if 

C  Assign  SP  scattering  matrix 

SP  (1, 1)  =  tSP  (1,1) 

SP  (2,1)  =  tSP  (2, 1) 
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SP  (1,2) 

= 

tSP 

(1, 

2) 

SP  (2,2) 

= 

tSP 

(2, 

2) 

Assign 

SS 

scatte 

ring  matrix 

SS  (1, 1) 

= 

tss 

(1, 

1) 

SS (2,  1) 

= 

tss 

(2, 

1) 

SS (1, 2) 

= 

tss 

(1, 

2) 

SS (2, 2) 

= 

tss 

(2, 

2) 

return 

end 

1 

end 

of 

subroutine 

subroutine  READ LAYER ( layer f 1 , NMESH, wnum, eOUT, ecode ) 

C 

C  reads  layer  parameters  from  data  file 

C 

implicit  none 
real*8  pi 

parameter  (pi=3 .1415 92 6535897932384 6d0) 
integer  LNR  !  max  number  of  resonances 
parameter  (LNR  =  100) 

integer  INTYPE  !  layer  input  type  wnum(I) 
integer  NMESH  !  number  of  frequencies 
integer  NRES  !  number  of  resonances 

real*8  wnumN (LNR) , sxN (LNR) , syN (LNR) , szN (LNR) , gN (LNR)  !  resonances 
real*8  EINF,CV0L  !  dielectric  constant  at  infinity  and  cell  volume 
real*8  NR,NI  !  index  of  refraction  (NR,NI) 
real*8  ER,  El  !  permittivity  E  =  (ER,  El) 

real*8  ERT1,  ERT2,  EIT1,  EIT2  !  permittivity  from  input  table  wnum(I) 

real*8  ERTD,EITD  !  permittivity  difference  in  input  table 

real*8  wnumTl ,  wnumT2  !  wave  numbers  from  input  table 

real*8  wnumTD  !  wave  numbers  difference  in  input  table 

integer  NTAB  !  number  of  entries  in  data  table 

integer  NTABERR  !  error  flag  for  mesh/data  table  mismatch 

real*8  wnum(*)  !  wavenumber 

integer  IL,I,J,JS  !  counters 

COMPLEX* 1 6  eOUT ( * ) 

character*79  layerfl  !  layer  data  file 
integer  ecode  !  errorcode 

COMPLEX* 1 6  dcmplx 
intrinsic  dcmplx 

ecode  =  0 

open (unit  =  10 ,  file=layerfl , status= ' old iostat=ecode )  !  layer  data 

if  (ecode  .eq.  0)  then 

read ( 10 , * )  INTYPE  !  ambient  region  input  type 
if  (INTYPE  .eq.  0)  then 

read  (10,*)  NR,  NI  !  index  of  refraction 
call  N2E (NR, NI, ER, El) 
do  1=1, NMESH 

eOUT  (I)  =  dcmplx (ER, El) 
end  do 

elseif  (INTYPE  .eq.  1)  then  !  dielectric  constant 
read (10, *)  ER,  El 
do  1=1, NMESH 

eOUT (I)  =  dcmplx (ER, El) 
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resonances 


end  do 

elseif  (INTYPE  .eq.  2)  then  ! 
read (10, *)  EINF,  CVOL 
read (10, *)  NRES 
print  *,  ' NRES= ' , NRES 

do  J=1 , NRES 

read ( 10 , * )  wnumN (J)  ,sxN (J) , syN ( J) , szN (J) ,gN(J) 
end  do 

call  SUMOFRES (ER,  El, EINF, CVOL, wnum, 

&  wnumN, sxN, syN, szN, gN, NRES) 

do  1=1 , NMESH 

eOUT(I)  =  dcmplx (ER, El) 
end  do 

elseif  (INTYPE  .eq.  3)  then  !  table 
do  1=1, NMESH 

if  (ecode  . ge .  0)  then 
NTABERR  =  1 
close  ( 10 ) 

open (unit =10 , f ile=layerf 1 , status='old' , iostat=ecode ) 
if  (ecode  .eq.  0)  then 

read (10,*)  INTYPE  !  input  type 

read (10,*)  NTAB  !  number  of  entry  lines 

JS  =  2  !  read  first  two  lines 

read (10, *)  wnumTl,  ERT 1 ,  EIT1 

read (10, *)  wnumT2,  ERT2 ,  EIT2 

do  J=3 , NTAB 

if  ((wnum (I)  .It.  wnumTl)  .or. 

&  (wnum (I)  . gt .  wnumT2) )  then 

wnumTl  =  wnumT2 
ERT 1  =  ERT2 
EIT1  =  EIT2 

JS=  JS  +  1  !  increase  line  number  by  1 

read ( 10 , * )  wnumT2,  ERT2 ,  EIT2 
end  if 
end  do 

if  (.not.  (((wnum (I)  .It.  wnumTl)  .or. 

&  (wnum (I)  . gt .  wnumT2) ) ) )  then 

wnumTD  =  wnumT2  -  wnumTl 
ERTD  =  ERT2  -  ERT 1 

EITD  =  EIT2  -  EIT1 

if  (wnumTD  . gt .  Id-300)  then 

ER  =  ERT1  +  ERTD* (wnum ( I )  -  wnumTl) /wnumTD 

El  =  EIT1  +  EITD* (wnum ( I )  -  wnumTl) /wnumTD 

NTABERR  =  0 
else 

ER  =  ERT1 

El  =  EIT1 

NTABERR  =  2 
end  if 
end  if 

eOUT(I)  =  dcmplx (ER, El ) 
if  (NTABERR  .eq.  1)  then 
ecode=2000 

print  *,  'INPUT  ERROR:  mesh/table  mismatch' 
end  if 

end  if  !  ecode  =  0 
end  if  !  ecode  >=  0 
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end  do 

end  if  !  INTYPE 
end  if  !  ecode  =  0  on  input 
close  (10)  !  close  layer  data 

return 

end  !  end  of  subroutine 


Brief  Tutorial 

1.  Setup  the  geometry  and  composition  of  layers  in  input  file  (MLAYER.INP). 

Example: 

0  1 

50.0  150.0d0  1001  1 

40.0d0  90.0d0  101  1 

2 

- - - ambient  layer 

data/ AIR 

— . - - - layer  1 

data/HMX.  table 
l.OdO 

- layer  2 

data/GOLD 

54.5d-l 

- - - layer  3 

data/ AIR 
l.OOdO 

The  first  line  is  for  units  (0  -  CGS;  1  -  SI)  and  scan  type  (0  -  angle  scan;  1  - 
wavenumber  scan). 

The  second  and  the  third  lines  contain  the  range  of  scan,  the  number  of  mesh 
points  and  the  bin  size:  the  second  line  is  for  the  wavenumber  scan  and  the  third 
line  is  for  the  angle  scan. 

The  fourth  line  contains  the  number  of  layers  (excluding  ambient  layer  and 
substrate). 

The  fifth  line  is  a  divider. 

The  sixth  line  is  for  the  path  and  the  file  name  of  the  data  file  of  the  ambient  layer. 
The  sixth  line  is  a  divider. 

The  seventh  line  is  for  the  path  and  the  file  name  of  the  data  file  of  the  first  layer. 
The  eighth  line  contains  the  thickness  of  the  first  layer. 

The  last  sequence  of  lines  (from  the  sixth  to  the  eighth  is  repeated  as  many  time  as 
many  layers  are  declared  on  the  fourth  line  plus  one  extra  set  for  the  substrate 
(layer  3  in  our  example). 

2.  Place  data  files  for  each  layer  at  the  location  specified  in  MLAYER.INP  file. 
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Examples  of  data  files: 

AIR  file 

0 

l.OdO  O.OdO 

The  first  line  is  for  input  type: 

0  -  real  and  imaginary  parts  of  the  refraction  index 

2  -  table  of  resonances  (frequency,  oscillator  strength,  broadening) 

3  -  real  and  imaginary  parts  of  the  permittivity  in  the  form  of  data  table 

HMX.table  file 


3 

201 

50.0000000000000 

50.4975124262273 

50.9950248524547 

51.4925372786820 

51.9900497049093 

52.4875622242689 

52.9850745573640 

53.4825868904591 

53.9800994098186 

54.4776119291782 


4.25905358483179 

4.28289690235066 

4.30819579494899 

4.33508465558413 

4.36371943716915 

4.39427844001136 

4.42696385313709 

4.46200369947143 

4.49965349284311 

4.54019690996128 


0.124136870938063 

0.131953983279890 

0.140717478521852 

0.150560832568167 

0.161650246864707 

0.174192141418067 

0.188442912146539 

0.204721994922112 

0.223429373616938 

0.245069178213272 


146.019899845123 

146.517413854599 

147.014927864075 

147.512435913086 

148.009949922562 

148.507463932037 

149.004977941513 

149.502485990525 


2.33534257788569 

2.34274505049977 

2.34997051451085 

2.35702540917843 

2.36391609417684 

2.37064836425402 

2.37722780923246 

2.38365967068105 


3.796482250904258E-002 
3.704172732575566E-002 
3.615343225924952E-002 
3.529812379956010E-002 
3.447408 11 4083763E-002 
3.367972744632545E-002 
3.291357929336934E-002 
3.217425750164293E-002 


The  first  line  is  for  input  type  (3  -  data  table) 

The  second  line  contains  the  number  of  data  points  in  the  table. 

The  table  contains 
wavenumber  (first  column), 

permittivity:  real  (second  column)  and  imaginary  (second  column)  parts. 
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HMX  file 


2 

2.9d0  l.OdO 
8 


46.5d0 

O.OOdO 

0.08d0 

O.OOdO 

5.00d0 

62.9d0 

O.OOdO 

61.54d0 

O.OOdO 

5.00d0 

76.3d0 

10.36d0 

O.OOdO 

5.22d0 

5.00d0 

84.2d0 

30.7  ldO 

O.OOdO 

29.97d0 

5.00d0 

87.9d0 

O.OOdO 

79.15d0 

O.OOdO 

5.00d0 

97.8d0 

O.OOdO 

29.38d0 

O.OOdO 

5.00d0 

100. 3d0 

4.75d0 

O.OOdO 

204.88d0 

5.00d0 

116.0d0 

57.01d0 

O.OOdO 

68.39d0 

5.00d0 

The  first  line  is  for  input  type  (2  -  table  of  resonances). 

The  second  column  contains  high-frequency  permittivity  and  the  volume  of  the 
unit  cell. 

The  third  line  contains  the  number  of  data  points  in  the  table. 

The  table  contains 
wavenumber  (first  column), 

oscillator  strength:  x  (second  column),  y  (third  column),  and  z  (fourth  column) 
components, 

broadening  of  each  resonance 


Appendix  2 

EM  Response  Calculations  using  DFT 

Presented  in  this  section  is  a  brief  tutorial  that  describes  the  sequences  of 
procedures  for  calculation  of  response  spectra  using  density  functional  theory  (DFT).  The 
simulation  framework  presented  here  adopts  the  DFT  software  NRLMOL.  In  that 
NRLMOL  is  a  general  purpose  DFT  model  for  application  to  many  different  types  of 
analysis,  this  tutoring  is  structured  to  put  emphasis  on  those  aspects  and  associated 
procedures  of  NRLMOL  that  are  important  for  the  numerical-simulation  framework 
presented  here. 

1.  fill  in  <name>. CLUSTER  file 
-number  of  atoms 

-positions  of  atoms  in  atomic  units/atomic  numbers 

2.  run: 

000_prepare_relax  <name> 
creates  directory  <name>. RELAX 
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010_submit_relax  <name> 
submits  relax 

check  GEOCNVRG  for  'CONVERGE  TRUE1  message 
020_copy_into_FREQ_subdir  <name> 
creates  directory  <name>.FREQ 
030_generate_with_specsym  <name> 

040_prepare_first  <name> 

050_submit_first  <name> 
check  for  conversion: 

'grep  SELF  nitro.specsym_l.stdouf  should  return  'SELF-CONSISTENCY  REACHED' 
060_prepare_all_the_rest 
07 0_submit_all_the_rest 
080_c  ombine_FRC  OUT 
090_final_spec  sym 
see  specsym.out  and  infred.spc  files 
1 00_prepare_frozen_phonon 
1 10_submit_frozen_phonon 
1 20_c  ombine_FRC  OUT_frozen_phonon 
130_specsym_frozen_phonon 
see  specsym.out. FPH 


27 


